Finite Temperature Critical Behavior of Mutual Information 
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We study mutual information for Renyi entropy of arbitrary index n, in interacting quantum 
systems at finite-temperature critical points, using high-temperature expansion, quantum Monte 
Carlo simulations and scaling theory. We find that for n > 1, the critical behavior is manifest 
at two temperatures Tc and nTc- For the XXZ model with Ising anisotropy, the coefficient of the 
area-law has a tint singularity, whereas the subleading correction from corners has a logarithmic 
divergence, with a coefficient related to the exact results of Cardy and Peschel. For T < nTc there 
is a constant term associated with broken symmetries that jumps at both Tc and nTc, which can be 
understood in terms of a scaling function analogous to the boundary entropy of Affieck and Ludwig. 
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The numerical study of entanglement in quantum sys- 
tems, through the entanglement entropy (EE) at zero 
temperature or mutual information (MI) at non-zero 
temperature, promises to be a new approach to quantify- 
ing properties of quantum phases that cannot be detected 
using traditional measures based on two-point correla- 
tion functions. It has already been used in one dimen- 
sional (ID) systems to identify the central charge in 
two dimensional (2D) systems to test the area law in the 
Heisenberg model [3|, Ij] , and to identify a topologically 
ordered spin liquid phase in a 2D spin model [5| . 

In ID, gapless systems described by conformal field 
theory show logarithmic violations of the area law 
However in 2D, the presence of an area law for a system 
such as the Heisenberg model implies that the existence 
of gapless modes does not necessarily lead to such a vio- 
lation. Similar area-law behavior is also observed in gap- 
less 2D bosonic theories while non-interacting fermions 
show a logarithmic violation presumably reflecting 
the infinite number of gapless modes associated with the 
fermi surface. The question of precisely which interacting 
many-body models have enough entanglement to violate 
the area law is important both for identifying new phases 
and for developing novel computational tools. 

Even with an area law, subleading corrections to the 
entanglement entropy, such as those associated with cor- 
ners, can show logarithmic divergence at quantum criti- 
cal points 0, 3 ■ While entanglement entropy at T = 
remains a key focus of current research, MI at non-zero 
temperature (which reduces to EE at T = 0) can also 
show universal critical behavior and has been a subject 
of both theoretical and computational @ studies. 

Prom a computational point of view, there is a clear 
need for new methods capable of studying EE or MI for 
large-scale quantum systems in D > 1. In this paper, we 
develop a High Temperature Expansion (HTE) method 
for calculating MI for Renyi entropy of arbitrary index 
n for lattice models in the thermodynamic limit. Our 
work represents a new direction in the use of series ex- 



pansions to study boundary phenomena in critical sys- 
tems, enabling one to calculate corner exponents such as 
Cardy-Peschel exponents in 2D systems [10| • 

In the following, we combine HTE with quantum 
Monte Carlo (QMC) simulations and a scaling theory 
to obtain the critical behavior of MI for a 2D spin-1/2 
XXZ model, H = E{^,J)iSfS^ + Sf + A5f5|), with 
Ising anisotropy A = 4. Quite generally, we find that 
for n > 1, the critical behavior manifests itself at two 
different temperatures Tc and nTc- Since this model is in 
the universality class of the 2D Ising model, the singu- 
larity of the area-law term is known to be t Int, where 
the reduced temperature t is |T — Tc\ or \T — nTc\, and 
the logarithmic divergence of the subleading corner terms 
can be related to the work of Cardy and Peschel • We 
also find that spontaneously broken symmetries lead to 
a constant term in the MI that jumps at nTc and Tc, de- 
scribed by a scaling function analogous to the boundary 
entropy of Affleck and Ludwig \vi\ . 

Replica Calculation of MI — Consider a system divided 
into two regions A and B, where pA is the reduced density 
matrix on A. The Renyi entropies are deflned as 

Sn{pA) = -^ln[Tiipl)]. (1) 

1 — n 

The von Neumann entropy 5*1 is defined by the limit n — s- 
1. The advantage of the n > 1 Renyi entropies is that 
they can be calculated by a "replica method" for integer 
n pj, where for a given inverse temperature /3, one must 
evaluate a partition function Z[A, n, T] corresponding to 
a path integral on a system with modified space-time 
topology. In region A, the system is periodic with period 
n/3, while in region B there are n distinct sheets, each 
periodic with period /?. Normalizing correctly, one has, 

Z[A,n,T] 
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where Z[T] denotes the partition function at temperature 
T. This replica method was used in for QMC simula- 
tions of 5*2. In this paper, we perform similar simulations 
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FIG. 1: In the HTE (left), we consider partitions of the 
infinite square-plane where region A could be a half-plane 
such as a U fe or a quadrant such as b or c. In QMC (right) the 
A''-site real space lattice is a torus with periodic boundaries, 
or a cylinder with the dashed boundaries open. In the torus, 
region A can be e ( "square" ) or e U / ( "strip" ) , while for the 
cylindrical we use a strip (e U /) region A. 



(for 71 = 2 and higher) and also develop a HTE method to 
calculate the partition function with this modified topol- 
ogy in powers of /3. We are interested in determining the 
MI between region A and its complement B, defined as 



Sn{PA) + Snips) - Snip). 



(3) 



One important feature of this calculation is that if the 
given Hamiltonian has a critical point at temperature 
Tc, then the partition function Z[A,n,T] shows critical 
behavior at T = nTc because the path integral is periodic 
with period n(3 in region A. If we consider a semi-infinite 
region A (for example, dividing a 2D plane into two half- 
planes) then Z[A,n,T] will be non-analytic at T = nTc- 
In contrast, the von Neumann MI, /i, should not show 
critical behavior at temperatures other than Tc- 

High Temperature Series and MI — We develop a HTE 
for the MI in powers of /3 = 1/T. One important simplifi- 
cation of MI is that all terms proportional to the volume 
oi A,B cancel out and we are left with only terms lo- 
calized near the boundary of A, B. The reason for this 
is that a given bulk term in region A appears once at 
temperature T/n in Z[A,n,T] but also appears once at 
temperature T/n in Z\T/ri\ - these cancel out. Similarly, 
this term appears n times in Z\B,n,T] but also appears 
n times in nln(Z[r]). 

The HTE is calculated by a linked cluster method 
12l |l3| . We imagine that the infinite system is divided 
into subregions A and B either by a single straight line 
running parallel to one of the axes, or by two perpendic- 
ular lines that meet at a point (Fig. 1). The line contri- 
bution to the MI is obtained by considering region A to 
be the half-plane aUb. To obtain the corner contribution 
we consider four separate partitions of the square lattice: 
the region A can be (i) the quadrant b (ii) the quadrant c 
(iii) the half-plane a U 6 or (iv) the half plane a U c. If we 
add MI from the first two partitions and subtract those 
from the next two, all line contributions cancel. The dif- 
ference defines two times the contribution from a single 
corner. More generally, we express /„ as 



TABLE I: High temperature series coefficients for the line 
and corner terms for Renyi MI for n = 2. 
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1.125 
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0.375 
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6.32421875 


-2.765625 


5 


-5.109375 


0.46875 


6 


64.02701823 


-27.11848958 


7 


15.59501953 


-0.3969401042 


8 


1079.586016 


-584.0700043 


9 


97.15596924 


-63.38234592 


10 


12847.34193 


-8700.183385 


11 


-1079.890682 


94.58389488 



where a„,6n,dn depend on /3, L is the length of the 
boundary. Tic is number of corners, and o?„ is a con- 
stant term, associated with symmetry breaking, to be 
explained later. 

Before division by the factor (1 — n), the coefficient of 
/3™ is a polynomial in ?i of order m, which vanishes at 
71 = and n ~ \. Thus dividing by 1 — 7i and taking the 
limit 71 — > 1 is simple and reduces the final coefficient to 
a polynomial of order 771 — 1. The complete expression 



for the line term to /3 is: 
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+71(71^ +71—1) 



C4 
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Here, A2 

v2 



2 + A2, A3 = -6A, A4 = (2 + A2)2 + 4(1 + 
2A^), B4 = A^- 4A2 and C4 = 2 + A". In addition, we 
have calculated both the line and corner contribution for 
the second Renyi entropy up to order Let 
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a„(/3) •i + ?ic6„(/3) + d„(/3). 



(4) 



The coefficients 1^ and up to m = 11 for the second 
Renyi entropy for A = 4 are given in Table ID 

Comparison with Exact Numerics - We calculate the 
MI via exact diagonalization (ED), and Stochastic Series 
Expansion [ij] QMC using the rephca-trick. Eg. We 
extend the QMC algorithm outlined in Ref. [9| to allow 
calculations to arbitrary Renyi entropies by directly con- 
structing a simulation cell with 71 sheets [l6| . Geometries 
considered are illustrated in Fig. 1. The critical temper- 
ature of the model is best determined by studying the 
Binder ratios associated with the order parameter. We 
estimate 2.234 < Tc < 2.237, which gives /3c/2 in the 
range 0.2235 to 0.2238. 
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FIG. 2: (color online) A comparison of MI calculated with 
HTE, and QMC on 32 x 32 simulation cells on toroid and 
cylindrical geometries (see Fig. 1). Inset: MI plotted against 
L at /3 = 0.204, showing excellent fit to linear form. 



Fig. [21 shows the partial sums of the 11-th order series 
for the hnear terms (area-law) for n = 2 compared with 
QMC data. One can see that the results agree extremely 
well up to a /3 value of 0.21, at which point the QMC 
data shows a sharp rise. To study the critical behavior 
more closely we use Pade approximants. The 2D Ising 
universality class is special in that the correlation length 
exponent v = 1 and the boundary free energy has a tint 
singularity [l5|. Anticipating this, we take two deriva- 
tives of the series, and use Pade approximants biased to 
have a pole at the /3c/2 value obtained from the Binder 
ratios. Upon integration, these lead to a tint singularity. 
One such approximant is shown in Fig. [2] (thick dashed 
line). It captures the sharp rise in QMC data extremely 
well, confirming the tint behavior to high accuracy. 

The corner terms should have a logarithmic singularity. 
In fact, the series for db2/d/3 show good convergence for 
a simple pole implying that 62 goes as some constant x 
times Int. To get an accurate estimate for the coefficient 
X, we once again bias the critical temperature values. 
With the critical point biased at 0.223 the spread of Pade 
approximants leads to an estimate of a; = 0.0143 ±0.0013 
, where as biasing it at 0.224 leads to an estimate of 
0.0151 ± 0.0016. We can relate these coefficients to the 
exact results of Cardy and Peschel [13] • The internal 
angle for the corner is 7 = tt/2 for region A and 7 = 37r/2 
for region B. Together with c ~ 1/2 for the Ising model, 
Eq. 4 in Rcf . [lO| leads to a —j^lnL singularity at the 
critical point. This, using v = 1, translates in to an x 
value of ^ = 0.0138. Our results show that for exactly 
soluble 2D universality classes with known values of the 
central charge, the results of Cardy and Peschel can be 
used to obtain the coefficient x. 

Fig. [3] shows a comparison of the von Neumann MI, 
calculated by continuing the HTE to n = 1, with results 



FIG. 3: (color online) MI divided by boundary length from 
the fourth-order series expansion for: Ji, compared to ED 
on a 4 X 4 system with 3x3 region A; and I4,, compared 
to QMC for 10 X 10 and 20 x 20. A sharp change in slope 
in I4/L occurs at l/4rc, where Tc = 2.24, obtained from the 
crossing of the fourth-order Binder cumulant for the staggered 
magnetization. Inset: the constant scaling term di extracted 
from a linear fit to the QMC data. 



obtained by exact diagonalization (ED) on a 4 x 4 sys- 
tem. In this case, we have multiplied the series by the 
length of the boundary separating regions A and B. The 
agreement is excellent up to /? « 0.25, which confirms the 
validity of both calculations and shows that finite size ef- 
fects are small at smaller /3 values. The von Neumann 
entropy series should be convergent down to Tc- Fig. [3] 
also compares HTE and QMC simulation results for I4, 
which further confirms that for /„ the higher tempera- 
ture singularity moves to nTc- The inset illustrates the 
constant scaling term c?4(/3) extracted from QMC data 
taken on 10 x 10 and 20 x 20 toroidal simulation cells 
with strip regions A (e U / in Fig. 1). At low temper- 
atures, di{l3) approaches the value ln(2) predicted from 
our scaling theory. For temperatures between Tc and 
nTc, theory predicts that dn = — ln(2)/(n— 1), discussed 
below, which is visible as a plateau in the QMC data. 

Constant Terms in the MI Due to Symmetry 
Breaking — We now consider the MI between region A 
and B away from criticality, in the limit of large sys- 
tem size. In addition to the line and corner terms, sym- 
metry breaking can lead to additional constant terms 
dn- First, consider the case of T < Tc, where the 
Ising symmetry is broken in all regions. The break- 
ing of the symmetry means that the partition functions 
Z[T],Z[T/n],Z[A, T, n],Z[B, T, n] aU have a multiplica- 
tive factor of 2 in addition to the volume, line, and cor- 
ner terms. The volume terms still cancel, and the line 
and corner terms still contribute according to Eq. ([4]), 
but the factors of 2 increase the MI by ln(2). Simi- 
larly, for Tc < T < nTc, the partition function Z(T) 
has no additional factors of 2 but Z{A, n, T), Z{B, n, T) 
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and Z{T/n) do, giving rise to a constant term in MI of 
dn = — ln(2)/(n — 1). These are verified by tfie plateau 
in the QMC data in the inset of Fig. [31 These resuhs 
are modified strongly by finite size effects due to the vol- 
ume terms which go generically as exp(— L/^), but can 
be substantial when ^ is comparable to or larger than L. 
These form part of the scaling theory, which we develop 
next. 

Scaling Theory Near nTc — Near T = nTc we can use 
scaling theory to describe the singular behavior of the 
Renyi entropy. The sheets of the system with period (3 
are not critical, while the region with period n(3 is in a 
critical scaling regime. Consider first the case that T > 
nTc and L ^ ^. Then, we can calculate the MI by using a 
scaling ansatz for the free energy of a critical theory with 
a boundary, which implies that the singular terms in the 
MI equal ci{L/^)+C2nc In(^) for some universal constants 
ci , C2 . The ci term represents the fact that the singular 
terms in the MI are due to degrees of freedom at length 
scale ^ and there are L/^ such terms. For T < nTc and 
L ^ ^, there is the additional — ln(2)/(n — 1) described 
above, but the singular MI behaves again as c[{L/^) + 
C2ncln(^). For L ^ ^, finite size scaling implies that 
the MI is equal to F{Lt^) + ric ln(min(L, ^)) plus smooth 
terms (such smooth terms multiplying L or nc) where 
F{x) cix as X +00, F{x) —?■ c[\x\ - ln(2)/(n - 1) 
asx— cx). Ata; = 0, F(x) equals n times the Affleck- 
Ludwig boundary entropy 

The 2D Ising universality class with = 1 is special 
and in this case the line term has a multiplicative log 
correction as verified in our series analysis. The sublead- 
ing corner term is predicted to diverge logarithmically, 
in agreement with the series calculation. The negative 
jump in the additive constant term together with an in- 
creasing line term leads to an approximate crossing of 
1,1 /L for different system sizes near nTc as seen in the 
QMC data in Fig. El 

Scaling Theory at Tc — At T near Tg, we can again de- 
velop a scaling theory. In contrast to the case of T « nTc, 
the region with period n/? is now in the ordered phase, 
and the n sheets with period f3 display critical scaling of a 
theory with a boundary magnetization (since the region 
with period n/S is ordered). For L ^ ^, the singular terms 
in the MI again behave as C3{L/£^)+C4nc ln(^)— ln(2) / {n— 
1) or c'^{L/^) + C4ncln(^) -|- ln(2) depending on whether 
T > Tc OT T < Tc with universal constants 03,03,04. 
The change in sign in the constant from — ln(2)/(n — 1) 
to ln(2) leads to a crossing of /„/£ for different system 
sizes at Tc (with corrections from line and corner terms 
which shift the crossings at finite L to larger T). There 
is again a multiplicative log correction in the Ising case. 
This critical point corresponds to the case of a bound- 
ary magnetic field, while the T = nTc critical point cor- 
responds to the case of free boundary conditions - but 
both produce a log correction. 

Discussion — We have developed computational meth- 



ods and scaling theory to study Renyi mutual informa- 
tion /„ in interacting quantum systems. Away from crit- 
ical points the MI consists of line terms (area- law), cor- 
ner terms, and constant terms coming from symmetry- 
breaking. At the critical points the line terms develop 
a singularity which vanishes as 1/^, and thus have a 
critical exponent v. In the special case of the 2D Ising 
universality class with v = 1, there are multiplicative 
log terms. The subleading corner terms show a log di- 
vergence, whose coefficient can be related to the central 
charge using the results of Cardy and Peschel [l^l- We 
also find that the constant terms jump discontinuously at 
the transitions and can be described by a scaling function 
that is analogous to the boundary entropy of Afficck and 
Ludwig 
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We have extended our previous QMC algorithm for 
I2 to calculate arbitrary by using a multi-sheeted 
space-time simulation cell, and confirmed the main re- 
sults of the scaling theory. QMC methods are able to 
access all temperature regions, allowing one to obtain 
the bulk terms due to symmetry breaking. HTE can 
separately obtain the line terms and subdominant corner 
terms. Since the HTE is immune to the sign problem, it 
should be a general tool for calculating MI in arbitrary 
interacting quantum systems such as frustrated spin or 
fermionic models in the future. 
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